\documentclass[10pt]{article}
\usepackage{amsmath}
\usepackage{bm}
\usepackage{bbm}
\usepackage{mathrsfs}
\usepackage{graphicx}
\usepackage{wrapfig}
\usepackage{subcaption}
\usepackage{epsfig}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{wrapfig}
\usepackage{graphicx}
\usepackage{psfrag}
\newcommand{\sun}{\ensuremath{\odot}} % sun symbol is \sun
\let\vaccent=\v % rename builtin command \v{} to \vaccent{}
\renewcommand{\v}[1]{\ensuremath{\mathbf{#1}}} % for vectors
\newcommand{\gv}[1]{\ensuremath{\mbox{\boldmath$ #1 $}}} 
\newcommand{\grad}[1]{\gv{\nabla} #1}
%\renewcommand{\baselinestretch}{1.2}
%\jot 5mm
%\graphicspath{{./figures/}}
%text dimensions
%\textwidth 6.5 in
%\oddsidemargin .2 in
%\topmargin -0.2 in
%\textheight 8.5 in
%\headheight 0.2in
%\overfullrule = 0pt
%\pagestyle{plain}
%\def\newpar{\par\vskip 0.5cm}
\begin{document}
%
%----------------------------------------------------------------------
%        Define symbols
%----------------------------------------------------------------------
%
\def\iso{\mathbbm{1}}
\def\half{{\textstyle{1 \over 2}}}
\def\third{{\textstyle{1 \over 3}}}
\def\fourth{{\textstyle{{1 \over 4}}}}
\def\twothird{{\textstyle {{2 \over 3}}}}
\def\ndim{{n_{\rm dim}}}
\def\nint{n_{\rm int}}
\def\lint{l_{\rm int}}
\def\nel{n_{\rm el}}
\def\nf{n_{\rm f}}
\def\DIV {\hbox{\af div}}
\def\GRAD{\hbox{\af Grad}}
\def\sym{\mathop{\rm sym}\nolimits}
\def\tr{\mathop{\rm tr}\nolimits}
\def\dev{\mathop{\rm dev}\nolimits}
\def\Dev{\mathop{\rm Dev}\nolimits}
\def\DEV{\mathop {\rm DEV}\nolimits}
\def\bfb {{\bi b}}
\def\Bnabla{\nabla}
\def\bG{{\bi G}}
\def\jmpdelu{{\lbrack\!\lbrack \Delta u\rbrack\!\rbrack}}
\def\jmpudot{{\lbrack\!\lbrack\dot u\rbrack\!\rbrack}}
\def\jmpu{{\lbrack\!\lbrack u\rbrack\!\rbrack}}
\def\jmphi{{\lbrack\!\lbrack\varphi\rbrack\!\rbrack}}
\def\ljmp{{\lbrack\!\lbrack}}
\def\rjmp{{\rbrack\!\rbrack}}
\def\sign{{\rm sign}}
\def\nn{{n+1}}
\def\na{{n+\vartheta}}
\def\nna{{n+(1-\vartheta)}}
\def\nt{{n+{1\over 2}}}
\def\nb{{n+\beta}}
\def\nbb{{n+(1-\beta)}}
%---------------------------------------------------------
%               Bold Face Math Characters:
%               All In Format: \B***** .
%---------------------------------------------------------
\def\bOne{\mbox{\boldmath$1$}}
\def\BGamma{\mbox{\boldmath$\Gamma$}}
\def\BDelta{\mbox{\boldmath$\Delta$}}
\def\BTheta{\mbox{\boldmath$\Theta$}}
\def\BLambda{\mbox{\boldmath$\Lambda$}}
\def\BXi{\mbox{\boldmath$\Xi$}}
\def\BPi{\mbox{\boldmath$\Pi$}}
\def\BSigma{\mbox{\boldmath$\Sigma$}}
\def\BUpsilon{\mbox{\boldmath$\Upsilon$}}
\def\BPhi{\mbox{\boldmath$\Phi$}}
\def\BPsi{\mbox{\boldmath$\Psi$}}
\def\BOmega{\mbox{\boldmath$\Omega$}}
\def\Balpha{\mbox{\boldmath$\alpha$}}
\def\Bbeta{\mbox{\boldmath$\beta$}}
\def\Bgamma{\mbox{\boldmath$\gamma$}}
\def\Bdelta{\mbox{\boldmath$\delta$}}
\def\Bepsilon{\mbox{\boldmath$\epsilon$}}
\def\Bzeta{\mbox{\boldmath$\zeta$}}
\def\Beta{\mbox{\boldmath$\eta$}}
\def\Btheta{\mbox{\boldmath$\theta$}}
\def\Biota{\mbox{\boldmath$\iota$}}
\def\Bkappa{\mbox{\boldmath$\kappa$}}
\def\Blambda{\mbox{\boldmath$\lambda$}}
\def\Bmu{\mbox{\boldmath$\mu$}}
\def\Bnu{\mbox{\boldmath$\nu$}}
\def\Bxi{\mbox{\boldmath$\xi$}}
\def\Bpi{\mbox{\boldmath$\pi$}}
\def\Brho{\mbox{\boldmath$\rho$}}
\def\Bsigma{\mbox{\boldmath$\sigma$}}
\def\Btau{\mbox{\boldmath$\tau$}}
\def\Bupsilon{\mbox{\boldmath$\upsilon$}}
\def\Bphi{\mbox{\boldmath$\phi$}}
\def\Bchi{\mbox{\boldmath$\chi$}}
\def\Bpsi{\mbox{\boldmath$\psi$}}
\def\Bomega{\mbox{\boldmath$\omega$}}
\def\Bvarepsilon{\mbox{\boldmath$\varepsilon$}}
\def\Bvartheta{\mbox{\boldmath$\vartheta$}}
\def\Bvarpi{\mbox{\boldmath$\varpi$}}
\def\Bvarrho{\mbox{\boldmath$\varrho$}}
\def\Bvarsigma{\mbox{\boldmath$\varsigma$}}
\def\Bvarphi{\mbox{\boldmath$\varphi$}}
\def\bone{\mathbf{1}}
\def\bzero{\mathbf{0}}
%---------------------------------------------------------
%               Bold Face Math Italic:
%               All In Format: \b* .
%---------------------------------------------------------
\def\bA{\mbox{\boldmath$ A$}}
\def\bB{\mbox{\boldmath$ B$}}
\def\bC{\mbox{\boldmath$ C$}}
\def\bD{\mbox{\boldmath$ D$}}
\def\bE{\mbox{\boldmath$ E$}}
\def\bF{\mbox{\boldmath$ F$}}
\def\bG{\mbox{\boldmath$ G$}}
\def\bH{\mbox{\boldmath$ H$}}
\def\bI{\mbox{\boldmath$ I$}}
\def\bJ{\mbox{\boldmath$ J$}}
\def\bK{\mbox{\boldmath$ K$}}
\def\bL{\mbox{\boldmath$ L$}}
\def\bM{\mbox{\boldmath$ M$}}
\def\bN{\mbox{\boldmath$ N$}}
\def\bO{\mbox{\boldmath$ O$}}
\def\bP{\mbox{\boldmath$ P$}}
\def\bQ{\mbox{\boldmath$ Q$}}
\def\bR{\mbox{\boldmath$ R$}}
\def\bS{\mbox{\boldmath$ S$}}
\def\bT{\mbox{\boldmath$ T$}}
\def\bU{\mbox{\boldmath$ U$}}
\def\bV{\mbox{\boldmath$ V$}}
\def\bW{\mbox{\boldmath$ W$}}
\def\bX{\mbox{\boldmath$ X$}}
\def\bY{\mbox{\boldmath$ Y$}}
\def\bZ{\mbox{\boldmath$ Z$}}
\def\ba{\mbox{\boldmath$ a$}}
\def\bb{\mbox{\boldmath$ b$}}
\def\bc{\mbox{\boldmath$ c$}}
\def\bd{\mbox{\boldmath$ d$}}
\def\be{\mbox{\boldmath$ e$}}
\def\bff{\mbox{\boldmath$ f$}}
\def\bg{\mbox{\boldmath$ g$}}
\def\bh{\mbox{\boldmath$ h$}}
\def\bi{\mbox{\boldmath$ i$}}
\def\bj{\mbox{\boldmath$ j$}}
\def\bk{\mbox{\boldmath$ k$}}
\def\bl{\mbox{\boldmath$ l$}}
\def\bm{\mbox{\boldmath$ m$}}
\def\bn{\mbox{\boldmath$ n$}}
\def\bo{\mbox{\boldmath$ o$}}
\def\bp{\mbox{\boldmath$ p$}}
\def\bq{\mbox{\boldmath$ q$}}
\def\br{\mbox{\boldmath$ r$}}
\def\bs{\mbox{\boldmath$ s$}}
\def\bt{\mbox{\boldmath$ t$}}
\def\bu{\mbox{\boldmath$ u$}}
\def\bv{\mbox{\boldmath$ v$}}
\def\bw{\mbox{\boldmath$ w$}}
\def\bx{\mbox{\boldmath$ x$}}
\def\by{\mbox{\boldmath$ y$}}
\def\bz{\mbox{\boldmath$ z$}}
%*********************************
%Start main paper
%*********************************
\centerline{\Large{\bf PRISMS-PF}}
\smallskip
\centerline{\Large{\bf Mechanics (Infinitesimal Strain)}}
\bigskip
Consider a elastic free energy expression of the form:
\begin{equation}
  \Pi(\varepsilon) = \int_{\Omega}  \frac{1}{2} \varepsilon:C:\varepsilon  ~dV - \int_{\partial \Omega}   u \cdot t  ~dS
\end{equation}
where $\varepsilon$ is the infinitesimal strain tensor, given by $\varepsilon = \frac{1}{2}(\grad u + \grad u^T)$,  $C_{ijkl}=\lambda \delta_{ij} \delta_{kl}+\mu ( \delta_{ik} \delta_{jl}+ \delta_{il} \delta_{jk} )$  is the fourth order elasticity tensor, ($\lambda$, $mu$) are the Lame parameters, $u$ is the displacement, and $t$ is the surface traction.

\section{Governing equation}
Considering variations on the displacement $u$ of the from $u+\alpha w$, we have
\begin{align}
\delta \Pi &=  \left. \frac{d}{d\alpha} \left( \int_{\Omega}    \frac{1}{2}\epsilon(u+\alpha w)  : C : \epsilon(u+\alpha w) ~dV  - \int_{\partial \Omega}   u \cdot t  ~dS \right) \right\vert_{\alpha=0}\\
&=  \int_{\Omega}   \grad w : C :  \epsilon  ~dV -  \int_{\partial \Omega}   w \cdot t~dS\\
&=  \int_{\Omega}   \grad w : \sigma  ~dV -  \int_{\partial \Omega}   w \cdot t  ~dS
\end{align}
where $\sigma = C :  \varepsilon$ is the stress tensor and $t=\sigma \cdot n$ is the surface traction.\\

The minimization of the variation, $\delta \Pi=0$, gives the weak formulation of the governing equation of mechanics:

\begin{align}
\int_{\Omega}   \grad w : \sigma  ~dV -  \int_{\partial \Omega}   w \cdot t  ~dS = 0
\end{align}

If surface tractions are zero: \\
\begin{align}
R &=  \int_{\Omega}   \grad w :  \sigma ~dV = 0 
\end{align}

\section{Terms for Input into PRISMS-PF}
In PRISMS-PF, two sets of terms are required for elliptic PDEs (such as this one), one for the left-hand side of the equation (LHS) and one for the right-hand side of the equation (RHS). We solve $R=0$ by casting this in a form that can be solved as a matrix inversion problem. This will involve a brief detour into the discretized form of the equation. First we derive an expression for the solution, given an initial guess, $u_0$:
\begin{gather}
0 = R(u) = R(u_0 + \Delta u)
\end{gather}
where $\Delta u = u - u_0$. Then, applying the discretization that $u = \sum_i w^i U^i$, we can write the following linearization:
\begin{equation}
\frac{\delta R(u)}{\delta u} \Delta U = -R(u_0) \label{matrix_eqn}
\end{equation}
The discretized form of this equation can be written as a matrix inversion problem. However, in PRISMS-PF, we only care about the product $\frac{\delta R(u)}{\delta u} \Delta U$. Taking the variational derivative of $R(u)$ yields:
\begin{align}
\frac{\delta R(u)}{\delta u} &= \frac{d}{d\alpha} \int_{\Omega}   \nabla w :C: \epsilon (u+\alpha w) ~dV  \bigg{|}_{\alpha=0} \\
&=  \int_{\Omega}   \nabla w :C: \frac{1}{2}\frac{d}{d\alpha}\left[ \nabla(u+\alpha w) + \nabla(u+\alpha w)^T \right] ~dV \bigg{|}_{\alpha=0}\\
&= \int_{\Omega}   \nabla w :C: \frac{d}{d\alpha}\nabla(u+\alpha w)  ~dV \bigg{|}_{\alpha=0} \quad (due ~to ~the ~symmetry ~of ~C) \\
&= \int_{\Omega}   \nabla w :C: \nabla w  ~dV 
\end{align}
In its discretized form $\frac{\delta R(u)}{\delta u} \Delta U$ is:
\begin{equation}
\frac{\delta R(u)}{\delta u} \Delta U = \sum_i \sum_j \int_{\Omega} \nabla N^i : C : \nabla N^j dV ~\Delta U^j
\end{equation}
Moving back to the non-discretized form yields:
\begin{equation}
\frac{\delta R(u)}{\delta u} \Delta U = \int_{\Omega} \nabla w : C : \nabla (\Delta u) dV
\end{equation}
Thus, the full equation relating $u_0$ and $\Delta u$ is:
\begin{equation}
\int_{\Omega} \nabla w : \underbrace{C : \nabla (\Delta u)}_{r_{ux}^{LHS}} dV = -\int_{\Omega}   \grad w : \underbrace{\sigma}_{r_{ux}} ~dV
\end{equation}
The above values of $r_{ux}^{LHS}$ and $r_{ux}$ are used to define the equation terms in the following input file: \\
$applications/mechanics/equations.cc$



\end{document} 
